#include "cem/boundary_condition_object.h"
#include "cem/boundary_condition_utils.h"

#include <libmesh/point.h>
#include <libmesh/vector_value.h>

namespace cem
{

using libMesh::DenseMatrix;
using libMesh::DenseVector;
using libMesh::Number;
using libMesh::Point;
using libMesh::Real;
using libMesh::RealVectorValue;

/*
 * 在面上施加 Neumann 边界条件：f · v_t 的边界积分项（复数）
 * 方案：通过幅值/方向 Provider 统一计算，避免积分循环内重复分支判断。
 */
void
NeumannBoundaryCondition::assemble_face(const BoundarySpec & spec,
                                        const ElectromagneticsSystem & system,
                                        const std::vector<std::vector<RealVectorValue>> & phif_v,
                                        const std::vector<Real> & jxwf,
                                        const std::vector<Point> & normals,
                                        const std::vector<Point> & qpoints,
                                        unsigned int n,
                                        DenseMatrix<Number> & /*ke*/,
                                        DenseVector<Number> & fe_local) const
{
    auto amplitude_provider =
        cem::boundary::make_amplitude_provider(spec, system.layers(), system.omega());
    auto direction_provider = cem::boundary::make_direction_provider(spec);

    const unsigned int n_qpf = jxwf.size();
    for (unsigned int qp = 0; qp < n_qpf; ++qp)
    {
        const RealVectorValue nrm = RealVectorValue(normals[qp]);
        RealVectorValue d_vec = direction_provider->compute_direction(nrm);
        const std::complex<double> g_qp = amplitude_provider->compute_at_point(qpoints[qp]);

        const std::complex<double> w = std::complex<double>(jxwf[qp], 0.0);
        for (unsigned int i = 0; i < n; ++i)
        {
            const RealVectorValue v_t = phif_v[i][qp] - (phif_v[i][qp] * nrm) * nrm;
            const double proj = (d_vec * v_t);
            fe_local(i) += Number(g_qp * (proj * w));
        }
    }
}

/*
 * 在面上施加 Robin 边界条件：alpha * (E_t, v_t) 的边界积分（复数）
 */
void
RobinBoundaryCondition::assemble_face(const BoundarySpec & spec,
                                      const ElectromagneticsSystem & /*system*/,
                                      const std::vector<std::vector<RealVectorValue>> & phif_v,
                                      const std::vector<Real> & jxwf,
                                      const std::vector<Point> & normals,
                                      const std::vector<Point> & /*qpoints*/,
                                      unsigned int n,
                                      DenseMatrix<Number> & ke,
                                      DenseVector<Number> & /*fe_local*/) const
{
    const std::complex<double> alpha(spec.value, 0.0);
    const unsigned int n_qpf = jxwf.size();
    for (unsigned int qp = 0; qp < n_qpf; ++qp)
    {
        const RealVectorValue nrm = RealVectorValue(normals[qp]);
        for (unsigned int i = 0; i < n; ++i)
        {
            const RealVectorValue phi_t_i = phif_v[i][qp] - (phif_v[i][qp] * nrm) * nrm;
            for (unsigned int j = 0; j < n; ++j)
            {
                const RealVectorValue phi_t_j = phif_v[j][qp] - (phif_v[j][qp] * nrm) * nrm;
                const double m_t = (phi_t_i * phi_t_j) * jxwf[qp];
                ke(i, j) += Number(alpha * m_t);
            }
        }
    }
}

/*
 * 在面上施加 Dirichlet（切向分量）边界条件：投影到切向空间后以最小二乘方式累积（复数）
 */
void
DirichletBoundaryCondition::accumulate_dirichlet_face(
    const BoundarySpec & spec,
    const ElectromagneticsSystem & system,
    const std::vector<std::vector<RealVectorValue>> & phif_v,
    const std::vector<Real> & jxwf,
    const std::vector<Point> & normals,
    const std::vector<Point> & qpoints,
    unsigned int n,
    const std::vector<libMesh::dof_id_type> & dofs,
    std::unordered_set<libMesh::dof_id_type> & dirichlet_rows_keep,
    std::unordered_map<libMesh::dof_id_type, std::complex<double>> & dirichlet_value_num,
    std::unordered_map<libMesh::dof_id_type, double> & dirichlet_value_den) const
{
    const double tol_trace = 1e-14;
    auto amplitude_provider =
        cem::boundary::make_amplitude_provider(spec, system.layers(), system.omega());
    auto direction_provider = cem::boundary::make_direction_provider(spec);
    const unsigned int n_qpf = jxwf.size();

    for (unsigned int i = 0; i < n; ++i)
    {
        double l2_tangent = 0.0;
        std::complex<double> proj_num(0.0, 0.0);
        for (unsigned int qp = 0; qp < n_qpf; ++qp)
        {
            const RealVectorValue nrm = RealVectorValue(normals[qp]);
            RealVectorValue d_vec = direction_provider->compute_direction(nrm);
            const std::complex<double> g_qp = amplitude_provider->compute_at_point(qpoints[qp]);

            const RealVectorValue phi_t = phif_v[i][qp] - (phif_v[i][qp] * nrm) * nrm;
            const double w = jxwf[qp];
            l2_tangent += (phi_t * phi_t) * w;
            proj_num += g_qp * ((d_vec * phi_t) * w);
        }
        if (l2_tangent > tol_trace)
        {
            const auto gid = dofs[i];
            dirichlet_rows_keep.insert(gid);
            dirichlet_value_den[gid] += l2_tangent;
            dirichlet_value_num[gid] += proj_num;
        }
    }
}

void
ImpedanceBoundaryCondition::assemble_face(
    const BoundarySpec & spec,
    const ElectromagneticsSystem & system,
    const std::vector<std::vector<libMesh::RealVectorValue>> & phif_v,
    const std::vector<libMesh::Real> & jxwf,
    const std::vector<libMesh::Point> & normals,
    const std::vector<libMesh::Point> & /*qpoints*/,
    unsigned int n,
    libMesh::DenseMatrix<libMesh::Number> & ke,
    libMesh::DenseVector<libMesh::Number> & /*fe_local*/) const
{
    auto amp = cem::boundary::make_amplitude_provider(spec, system.layers(), system.omega());
    std::complex<double> Z = amp->compute_scalar();
    if (std::abs(Z) < 1e-18)
        return;
    std::complex<double> Y = 1.0 / Z;

    const double omega = system.omega();
    const double Yr = std::real(Y);
    const double Yi = std::imag(Y);

    const unsigned nq = jxwf.size();
    for (unsigned qp = 0; qp < nq; ++qp)
    {
        const RealVectorValue nrm = RealVectorValue(normals[qp]);
        for (unsigned i = 0; i < n; ++i)
        {
            const auto phi_t_i = phif_v[i][qp] - (phif_v[i][qp] * nrm) * nrm;
            for (unsigned j = 0; j < n; ++j)
            {
                const auto phi_t_j = phif_v[j][qp] - (phif_v[j][qp] * nrm) * nrm;
                const double m = (phi_t_i * phi_t_j) * jxwf[qp];
                // 复数矩阵条目：A + i B，其中
                // A = -ω·Yi·M_t，B = ω·Yr·M_t
                ke(i, j) += Number((-omega * Yi) * m + std::complex<double>(0.0, omega * Yr) * m);
            }
        }
    }
}

// 工厂实现保持不变
std::unique_ptr<BoundaryConditionObject>
make_boundary_condition_object(BoundaryType type)
{
    switch (type)
    {
        case BoundaryType::Neumann:
            return std::make_unique<NeumannBoundaryCondition>();
        case BoundaryType::Robin:
            return std::make_unique<RobinBoundaryCondition>();
        case BoundaryType::Impedance:
            return std::make_unique<ImpedanceBoundaryCondition>();
        case BoundaryType::Dirichlet:
            return std::make_unique<DirichletBoundaryCondition>();
        default:
            return nullptr;
    }
}

} // namespace cem